function dx = dif_orbit_dynamics(t,x)

global miu;

rx = x(1);
ry = x(2);
rz = x(3);
Vx = x(4);
Vy = x(5);
Vz = x(6);
r= sqrt(rx^2 + ry^2 + rz^2);

drx = Vx;
dry = Vy;
drz = Vz;
dVx = -miu * rx / r^3;
dVy = -miu * ry / r^3;
dVz = -miu * rz / r^3;
dx = [drx; dry; drz; dVx; dVy; dVz];